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SUMMARY 


This paper discusses incompressible Navier-Stokes solution methods with an emphasis on 
the pseudo compressibility method. A steady-state flow solver based on the pseudocom- 
pressibility approach is then described. This flow-solver code has been used to analyze 
the internal flow in the Space Shuttle main engine hoi-gas manifold. Salient features as- 
sociated with this three-dimensional realistic flow simulation are discussed. Numerical 
solutions relevant to the current engine analysis and the redesign effort are discussed along 
with experimental results. This example demonstrates the potential of computational fluid 
dynamics as a design tool for aerospace applications. 


KEY WORDS Incompressible Flow, Navier-Stokes Equations, Finite Differences 

INTRODUCTION 

Incompressible flows are encountered in many realistic engineering problems related to 
hydrodynamics, and also to low-speed aerodynamics. However, to date, computational 
flow simulation has not been a critical element in resolving many of these problems. For 
instance, impellers, automobiles, submarines, chemical reactors, etc., have been designed 
reasonably well by empirical means. As technologies advance, the design of modern flow 
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devices lends to be more compact and highly efficient. Therefore, the approach relying on 
empiricism and simplified analysis becomes inadequate for resolving problems associated 
with those devices requiring advanced analysis. For example, in analyzing and redesigning 
the current Space Shuttle main engine (SSME) power head, computational simulations 
became an economical and time-saving supplement to experimental data. There are vast 
numbers of other real-world problems which demand accurate viscous, incompressible flow 
solutions, such as rocket-engine fuel flow, flow through an impeller (as in the turbo pump in 
the SSME), and biomedical flow (as in an artificial heart). There have been many computa- 
tional simulations reported in the past on these subjects. In response to these demands the 
present study attempts to develop an alternative method for simulating three-dimensional 
(3-D) viscous incompressible flows. This paper is a summary of our recent effort (NASA 
Ames, Rocketdyne, and NASA Marshall Space Flight Center) to develop and apply an 
incompressible Navier-Stokes solver with a special emphasis on SSME applications. 

A major problem involved in solving the incompressible Navier-Stokes equations comes 
from the lack of a pressure term in the continuity equation. In two-dimensional (2-D) 
problems, one can bypass this difficulty by using a stream function-vorticity formulation. 
In 3-D, a similar approach can be adopted using a vorticity-velocity formulation. However, 
the pressure term is removed at the expense of introducing three vorticity equations as 
well as requiring vorticity boundary conditions. In 2-D, computational efficiency is not a 
major problem^ however, in realistic 3-D applications, satisfying continuity in a reasonable 
amount of computing time becomes a primary issue. Naturally, computational efficiency is 
of primary" importance in addition to accuracy and robustness. For convenience anti flexi- 
bility, a primitive variable formulation in generalized curvilinear coordinates is chosen, and 
our discussion is limited to this formulation using a finite-difference approach. Derivation 
of equations and algorithmic details can be found in our earlier publications 1-3 and only 
equations relevant to the present discussion will be given in this paper. 

SOLUTION METHODS 


Unsteady, 3-D, incompressible flow" with constant density is governed by the following 
Navier-Stokes equations: 
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where t is the time, x, the Cartesian coordinates, u, the corresponding velocity components, 
p the pressure, and t (J the viscous-stress tensor. Here, all variables are nondimensionalized 
by a reference velocity and length scale. The viscous-stress tensor can be written as 


Tjj — 2 uSjj -R / j 

(3) 

1 du, Oil, 

(4) 

S ' J ~2 i S n + d r , ) 


Here, v is the kinematic viscosity, S,j is the strain-rate tensor, and R,j is the Reynolds 
stresses. Various levels of closure models for R, } are possible. In the present study, 
turbulence is simulated by an eddy viscosity model, using a constitutive equation of the 
following form: 

R;j = ^Rkkbjj - 2i' t Sjj (5) 

where, v t is the turbulent eddy viscosity. By including the normal stress, R a., in the 
pressure, u in equation (3) can be replaced by (u + rq), namely, 

Tjj = 2( v 4 v, ) S , j - ‘lv T S u ( 6 ) 

]n the remainder of this paper the total viscosity, i ' T ■, will be represented simply by r. 
Therefore, the incompressible Navier-Stokes equations are modified to allow variable vis- 
cosity in the present formulation. 

In this section, a few solution methods used in previous flow solver development work are 
briefly discussed. 

Method using Poisson’s equation for pressure 

Perhaps the oldest and most commonly used method of satisfying the continuity equation 
is with the use of a derived equation, namely, Poisson’s equation for pressure. 4 The usual 
computational procedure is to solve the pressure field such that continuity is satisfied at the 
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next time level. This procedure usually requires a relaxation scheme iterating on pressure 
until the divergence-free condition is satisfied within a specified numerical accuracy. To 
obtain a divergence free velocity field numerically, spacial differencing of the Laplacian 
operator in Poisson’s equation has to be consistent with the velocity differencing. Instead 
of a Laplacian operator, divergence of a gradient operator may be necessary. Even though 
this approach in general needs acceleration schemes to further enhance computational 
efficiency for 3-D problems, successful computations have been made. 5,6 

Fractiona l-step method 

The fractional-step method can be used for time-dependent computations of the incom- 
pressible Navier-Stokes equations. 7-14 Here, the time evolution is approximated by several 
steps. Various operator splitting can be adopted by treating the momentum equation as a 
combination of convection, pressure, and viscous terms. The common application of this 
method is done by two steps. The first step is to solve for an auxilary velocity field using 
the momentum equation in which the pressure-gradient term can be computed from the 
pressure in the previous time-step 10 , or can be excluded entirely. 11 In the second step, the 
pressure is computed which will map the auxiliary velocity onto a divergence-free velocity 
field. 

A generalized flow solver based on this approach using a staggered grid has been de- 
veloped in conjunction with the present work mainly for time-dependent computations. 
Details of the solver, which require lengthy description of various aspects such as the 
spatial-differencing scheme which maintains conservativeness, time integration and associ- 
ated boundary conditions, can be found in Reference 15. 


Pseudocompressibility method 


Recent advances in the state of the art in computational fluid dynamics (C'FD) have been 
made in conjunction with compressible flow computations. Therefore, it is of significant 
interest to be able to use some of these compressible flow algorithms. To do this, the arti- 
ficial compressibility method of C’horin 16 can be used. In this formulation, the continuity 
equation is modified by adding a time-derivative of the pressure term resulting in 


1 dp duj 

- — -f- — — 0 

d di d.v, 


( 7 ) 
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Together with the unsteady momentum equations, this forms a hyperbolic-parabolic type 
of time-dependent system of equations. Thus fast, implicit schemes developed for com- 
pressible flows, such as the approximate-factorization scheme by Beam and Warming 14 
or Briley and McDonald, 15 ’ 16 can be implemented. Various applications evolved from 
this concept have been reported for obtaining steady-state solutions. ” Merkle and 
Athavale 19 have successfully computed time-dependent flow's by subcycling each time step 
to satisfy the continuity equation. 

Strigburger, 20 and independently, Briley et ah, 21 interpreted the pseudocompressible for- 
mulation as a simplified form of low Mach number compressible equations. Then the low' 
Mach number approximation of the unsteady compressible flow equations becomes 


dp 

~di 


1 T du; dp 

T ( T7 7j * T u i ~n~~ = 0 

M dxj dxj 


( 8 ) 


where M is the Mach number. Therefore, the pseudocompressible formulation represented 
by equation (7) does not fully satisfy the time-dependent, compressible-flow formulation, 
and is intended mainly for steady-state computations. As shown above, unsteady-flow’ 
calculations may be performed by including the third term in the above equation. II 
approximate factorization is applied, where 1/(M 2 ) corresponds to j3, a large value of ,3 
will introduce factorization error. Therefore, the unfactored form would be preferable for 
time- dependent calculations. 


In an incompressible flow', a disturbance in the pressure causes waves which tiavcl with 
infinite speed. When pseudocompressibility is introduced, waves of finite speed result in 
which the magnitude of the speed depends upon the pseudocompressibility constant [3. In 
a true incompressible flow, the pressure field is affected instantaneously by a disturbance 
in the flow, but with pseudocompressibility, there will be a time lag between the flow 
disturbance and its effect on the pressure field. In viscous flow's, the behavior of the 
boundary layer is very sensitive to the streamwise pressure gradient, especially when the 
boundary layer is separated. If separation is present, a pressure wave traveling with finite 
speed will cause a change in the local pressure gradient which will affect the location of 
the flow separation. This change in separated flow will leed back to the pressuie field, 
possibly preventing convergence to a steady state. An extensive mathematical account on 
the pseudocompressibility approach is given by Temam." 


Among the various solution techniques, the method of of pseudocompressibility is chosen 
to gain efficiency and robustness in solving 3-D, real-world problems which require large 
numbers of grid points in curvilinear coordinates. In the present paper, the algorithm 
and its physical interpretation are described. Also, other practical considerations such 
as grid-induced error, boundary conditions, and a multiple-zone approach are presented. 
To improve the computational efficiency and accuracy, various acceleration schemes and 
modifications have been tried. The numerical scheme implemented in the present flow- 
solver (INS3D code) will be summarized next. 

NUMERICAL SCHEME 
Governing equations in generalized coordinates 

To perform calculations on 3-D, arbitrarily shaped geometries, the physical coordinates are 
t ransformed into general curvilinear coordinates by int roducing the following independent 
variables: 

T = f, (i = 

Applying the transformation to the governing equations (2) and (7) yields 

JLD+JL ( E._E vi) -_ = o (9) 

or d(i 

P' 

u 

V 

w . 

/«■'.■ + (fcMp-/*r 

uUj 4 - (£, ) x p 

vUi + (i,)y P 

wU,+U,),P -I 

and w-here J is the Jacobian of the transformation. The contravariant velocities, (with- 
out metric normalization), are defined as 



D _ 1 
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The viscous terms are given by 





o 1 



- 0- 


d_c. 

d.rb 


d , 

u 

du dv die 

+ ^ ix dr j + ^ y w J + ^ z 9 (j ) 


v ‘ 

V 

m w m 

■ 

i 





; 


When v is constant, the contribution of the second group of terms sums up to be zero 
when the velocity field is divergence- free. However, since in general, v varies in space and 
time, these terms have to be kept. For the flow with constant v in orthogonal coordinates, 
the above full viscous terms can be simplified further as follows: 




O' 

U 

V 

LIT. 


(116) 


Difference equations 

There are a number of different ways of defining variables in a grid system. Either standard 
cell node oriented or a staggered arrangement can be chosen. In Cartesian coordinates, a 
staggered grid has several favorable properties. 4,11 In generalized coordinates, these advan- 
tages become obscured because of the interpolation required. However, a fully conservative 
differencing scheme can be devised maintaining the convenience of a staggered arrange- 
ment in a Poisson solver. 12 Using any grid system, spatial differencing can be done either 
in finite-difference or finite- volume form. The finite- volume scheme usually behaves better 
near geometric singularities. In spacial differencing, both central differencing and upwind 
differencing have been tried. In this report, only the work using central differencing will 
be discussed. 


The numerical algorithm used to advance equation (9) in time is an implicit, noniterative, 
approximately factored, finite-difference scheme. 14 16 The time-differencing used by this 
scheme is generally known as the trapezoidal rule and is given by 


D n ] 1 - U u + 


Ar 

~2~ 


d by 


I <w 
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/ V dr 
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where the superscript n refers to the n 1h time step. 


The problem is to solve for D rt43 , and this is nonlinear in nature since £" 41 = E(D n *') 
is a nonlinear function of D" +1 . The following linearization procedure is used. A local 
Taylor expansion about 11 " yields 


E”+' = E? + - D n ) + 0(Ar 2 ) 

where .4, is the Jacobian matrix 


Ai = 


dEh 

dD 


(13) 


(14a) 


The Jacobian matrices can be represented by the following: 


4 - j 


rr 0 

I ii3 

L 2 j3 

L 3 }3 - 


Q -+■ L\ u 

L 2 « 

L 3 u 

Li 

L j r 

e + 

L 3 v 

-L 3 

L i w? 

L 2 w 

Q + L 3 «’- 


(146) 


where 

Q - L 0 + L \ « + L 2 v + L$w 

Lu — (£;)*> L] = (£t)x> -^2 = (£i)y> -^3 — {£i)z 
6 = (£,»/• «*() 

Substituting equation (13) into equation (9) results in the governing equation in delta form 


{ 


^ V 7 


n 4 1 


/ J" 41 

4 


where 


- r,) + ^(. 4 ; - r,>) + - r 3 )J j - 0 ") 

~ £>i ) n + bi){Eo - -Epi) 77 + ^'(^3 ■ £* p3 )" 

rr >)- 


r,/r 4! - £” 4 1 

h = At for trapezoidal 
= 2 ^t for Euler 

^ ~ finite - difference operator for , 


(is; 


o 


dt 
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At this point it should be noted that the notation of the form [^(j 4 - T)]D refers to 
%(AD) - $(TD) and not §f D - % D. 

Approximate Factorization 

The solution of equation (15) would involve a formidable matrix-inversion problem. With 
the use of an alternating direction implicit (ADI) type scheme, the problem could be 
reduced to the inversion of three matrices of small bandwidth, for which there exist some 
efficient solution algorithms. The particular ADI form used here is known as approximate 
factorization (AF). It is difficult to implement the AF scheme to equation (15) in its full 
matrix form. Noting that at steady state the left-hand side of equation (15) approaches 
zero, a simplified expression for the viscous term as shown in equation (lib) is used on the 
left-hand side. To maintain the accuracy of the solution, the entire viscous term is used 
on the right-hand side. Using this, the governing equation becomes 

LiL v L<(D n+ ' - D n ) = RHS (16) 

where 

L ( = / + ^J 7,+1 M^-'h) 

L v = /+ - 1 2 ) (17) 

i + ^> + \-U J-75) 

and RHS is the right-hand side of equation (15). When second-order central differencing is 
used, the solution to this problem becomes the inversion of three block tridiagonal matrices. 
The inversion problem is reduced to the three inversions 

(L n )±D =■- RHS 

(I<)A/> AD (18) 

(Z C )AZ) ,,+1 = A !) 

These inversions are carried out for all interior points, and the boundary conditions can 
be implemented explicitly. It is possible, however, to implement the boundary conditions 
implicitly. 
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The factorization has introduced the following second-order, cross-product term into the 
equation: 

h 2 (MiM* + M 2 M 3 + MsM 1 ] AD + o(h 3 ) (19) 

where 

v 4 i = A" - 71, A2 = A% - 72) ^3 = - 4 " - 73) h = 

To maintain the second order accuracy of the scheme, the added terms must be kept 
smaller than the original terms in the equations everywhere in the computational domain. 
This puts a restriction on the size of the pseudocompressibility parameter (3. The proper 
choice of /3 is discussed in our earlier reports. 1 ’ 2 In applying the approximate factorization 
scheme, it has been found that the stability of the scheme is dependent on the use of some 
higher-order smoothing terms. These are used to damp out higher-frequency oscillations 
which arise in the solution because of the use of second-order central differencing, and will 
be discussed in a following section. 

Numerical dissipation/smoothing 

Higher-order smoothing terms are required to make the present algorithm stable. These 
terms help to damp out the higher-order oscillations in the solution that are caused by 
the use of second-order central differencing. The smoothing term can be related to an 
upwind finite-difference approximation. An extensive discussion on numerical dissipation 
can be found in Reference 23. In this report, specifics relevant to the pseudocompressible 
formulation are discussed. 

Including these smoothing terms, equations (16) and (17) become 

- D n ) = RIIS of(15) - c £ .[(V<A<) 2 + (V.A,,) 2 + (V C A c ) 2 ]D n (16') 

where 

L ( = /+~J" +, » ( (.4r V ( A f 
i, = / + + -■>;,) + €, V,A, (17') 

Li --- I + y J " + '*< {A " - y,) + “ V (*< 
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Here, V and A represent forward and backward spacial-differencing operators, respectively. 
To preserve the tridiagonal nature of the system, only second-order smoothing can be used 
on the left-hand side of the equation, whereas fourth-order smoothing is used on the right- 
hand side. When the diagonal algorithm (described in a later section) is used, however, 
it is feasible to increase the bandwidth of the system to a pentadiagonal. This makes il 
possible to use fourth-order smoothing on the left-hand side of the equations also. W hen 
this combination of dissipation terms is used, the AF algorithm will be stable if t, and t c 
satisfy a certain relation. 23 

To study the nat ure of the numerical smoothing, one-dimensional form of dissipation terms 
are represented as follows: 

[1 - < : V ( A ( ] (p” +1 -p")= -MV{ A ()V (20) 

Suppose p is represented by the discrete Fourier expansion 

P = Y,mc iH ( 2 ') 

n 

where 


k — n = wave number 

A A(, 

77 = — A/2, ...0,1,. .(A/2 - 1) 

A = number of mesh poinis 

Substituting equation (22), equation (21) can be written as 

= ( 22 ) 

where 

k 9 — —2 + 2 cos(k) 

(k ! ) 2 — 6 — Scos(k) -f 2co.s(2A-) 

From this, the amplification factor can be defined as 

" = T"~ = — ir^Te] M) 
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To damp out the numerical fluctuations as time advances, the absolute value of the am- 
plification factor cr lias to be less than one for all possible frequencies, i.e. 

\<r\ < 1 

Noting that k 1 is always negative, this requirement leads to the following relation: 

c 6 < 2(1 - tils') (24 a) 

It can be shown that the above inequality is always satisfied if 

2c e < €j (246) 


The exact relation between these two coefficients can be determined only by a nonlinear 
stability analysis. In the present code, c, is taken to be three times larger than e € ,. From the 
expression given in equation (23), it is clear that if c? is too large, the rate of damping will 
be diminished. It may not be advantageous to take an excessively large value for over 
The choice of e t . depends on the Reynolds number and the grid spacings. However, as 
discussed later, large values of t € adversely affect the accuracy of the continuity equation. 
Therefore, the magnitude of e c is usually taken to be small. If grid sizes are fine enough 
to resolve the changes in the flow field, it can be taken as small as 10 -3 . 


There are two major sources of inaccuracy associated with the numerical dissipation terms; 
namely, (1) the numerical dissipation terms effectively change the Reynolds number of 
the flow, and (2) the explicit smoothing terms added to the continuity equation do not 
conserve mass. In particular, the explicit smoothing on the pressure can affect whether 
or not the solution converges to an incompressible flow field. Chang and Kwak 2 showed 
that the pseudo-pressure waves decay exponentially with time, and vanish as the solution 
converges. Thus the change in pressure with time approaches zero. When there is no 
explicit smoothing added to the equation, the divergence of the velocity field identically 
approaches zero. However, when explicit smoothing is included, as the change in pressure 
approaches zero, the divergence of velocity approaches 


du-i 

dx t 


MI) 

ft At 


[(v^) 2 + (v„d, ; ) 2 + (v^o 2 !/> 


(25) 
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where e f .(l) is the explicit smoothing parameter for the pressure. If t e is scaled by h, i.e., 
e ( = Ars e , equation (25) becomes independent of time step. This term can deteriorate 
the conservation of mass depending on the magnitude of j3 and the local pressure gradient. 
When the pressure gradients become substantial, as in the case when a region of sepa- 
ration is present, the smoothing term does not approach zero, and this contaminates the 
divergence of the velocity field. In this situation, mesh refinement usually does not help 
reduce the magnitude of this term. To allievate this problem, it may become necessary to 
decrease the smoothing coefficient in the continuity equation as the solution converges. 

Diagonal algorithm 

In a diagonal algorithm, a similarity transform is used to diagonalize the Jacobian matri- 
ces and uncouple the set of equations. The equations can then be solved by solving scalar 
tridiagonal matrices instead of solving block tridiagonal matrices. A similarity trans- 
form which symmetrizes and diagonalizes the matrices of the compressible gas dynamic 
equations has been used by Warming et al. 2 ^ and Turkel.” This method was used by Pul- 
liam and C’haussee 26 to produce a diagonal algorithm for the Euler equations. This can 
be applied to the compressible Navier-Stokes equations to obtain a considerable savings 
in computing time. 27 In this section, similarity transforms for the matrices used in the 
method of pseudocompressibility are presented. They are used in a diagonal algorithm 
which results in a substantial reduction in computer time. 

Similarity transformations exist which diagonalize the Jacobian matrices 

.4, = T,\,T-' (26) 

where A, is a diagonal matrix whose elements are the eigenvalues of the Jacobian matrices 
which is given by 

- Q 0 0 0 

- 0 Q 0 0 

Ai " 0 0 Q + C 0 

. 0 0 0 Q - C 

and where f'is the pseudospeed of sound, which is given by 

< ' = y (f? - Tq ) 2 + S(L] + L\ + L \ ) 
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The T, matrix is composed of the eigenvectors of the Jacobian matrix. For i=l (£-sweep), 
the first two eigenvectors are given by 
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For \—2 (7/-sweep), they are 
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For i = 3 ((.'-sweep), they are 


' o ■ 
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y>> 


k 

- Z v . 


- k - 


The determinant of T, is given by 

det(Ti) = 2 C 3 


(28) 


(29) 


(30) 


which remains bounded independent of the geometry. For more detail on the derivation of 
these matrices, see Reference 3- 


The implementation of the diagonal scheme involves replacing the Jacobian matrices in the 
implicit operators with the product of the similarity-transform matrices and the diagonal 
matrix as given in equation (26). The identity matrix in the implicit operators is replaced 
by the product of the similarity-transform mat rix and its inverse. A modification is made 
to the implicit viscous terms so that the transformation matrices may also be factored out 
of these terms. This implicitly adds an additional viscous dissipation term to the pressure. 
The transformation matrices are now factored out of the implicit operators to give 

C( = T ( [I 

c, = T,V+Y Jf M,-ii))T,' < 3I > 

£< = r<|7+^JMA<--)3)]r ( -‘ 
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where the implicit viscous terms are now given by 

7i = 

Since the transformation matrices are dependent on the metric quantities, factoring them 
outside the difference operators introduces an error. No modification has been made to 
the right-hand side of the equation; therefore, these linearization errors will not affect the 
steady-state solution, only the convergence path toward the solution. 

The implementation of this algorithm over the block algorithm will result in a substantial 
reduction in computational time per iteration because of the decrease in the number of 
operations performed. Additionally, considerably less memory is required to store the 
elements on the left-hand side. This additional memory was used to further vectorize 
the existing code as follows. Since the solution of a tridiagonal block or scalar matrix 
is recursive, it is not vectorizable for loops which use the current sweep direction as the 
inner do-loop index. However, if a large number of these matrices are passed into the 
inversion routines at once, then vectorization can take place in the ‘non-sweep’ direction. 
This diagonalized version is given as an option in the INS3D code. 

Boundary conditions 

An important part of any computer code is the proper implementation of boundary con- 
ditions. The code must be capable of handling the several different types of boundaries 
encountered in numerical simulations, which include solid-surface, in-flow and out-flow, 
and far-field boundaries. 

Solid surface 

At a solid-surface boundary, the usual no-slip condition is applied. In general the grid 
point adjacent to the surface will be sufficiently fine so that constant pressure normal to 
the surface in the viscous boundary layer can be assumed. For a (, = constant surface, this 
can be expressed as 



In-flow, out-flow and far-field conditions 


The in flow and the out-flow boundary conditions for an internal flow problem and far- 
field boundary conditions for an external flow problem can be handled in much the same 
way. The incoming flow for both problems can be set to some appropriate constant as 
dictated by the problem. For example, at the inlet to a pipe, the pressure can be set to 
a constant, and the velocity profile can be specified to be uniform. The conditions at the 
downstream, however, are the most difficult to provide. Simple upwind extrapolation is not 
well-posed. The best convergence rate is obtained if global mass is conserved. So to ensure 
the best results, the velocities and pressure are first updated using a second-order upwind 
extrapolation. Then these extrapolated velocities axe integrated over the exit, boundary 
to obtain the outlet mass flux. The extrapolated velocity components are weighted by the 
mass-flux ratio to conserve global mass, i.e., 





V * 


^ouf 


(33) 


If nothing further is done to update the boundary pressure, this can lead to discontinuities 
in the pressure because momentum is not being conserved. A method of weighting the 
pressure by a momentum correction is used to obtain a pressure which corresponds to the 
mass- weighted velocities 


P 


n+l 



1 

cl 


(trW r ) n+I 


(wWf] + ~(VC*VC) 


( (lw \ 71 f 1 / dw \ 


(34) 


where W is the contravariant velocity. In obtaining this formula, it has been assumed 
that the streamlines near the exit plane are nearly straight. Any appreciable deviation 
will cause a discontinuity in the pressure, and may lead to an instability. To avoid this, a 
momentum- weighted pressure was used. This was obtained by integrating the momentum- 
corrected pressure, p" +1 , and the extrapolated pressure, p", across the exit 

/; +1 = / p n+1 dd 

J tI H (.35) 

r p =f p* da 

J t* r. i t 

The final outlet pressure is then taken to be 

( 36 ) 
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Using these downstream boundary conditions, global conservation of mass and momentum 
are ensured, and the scheme will not introduce instabilities into the flow field. 

SSME POWER HEAD FLOW SIMULATION 
Validation of flow solver 

The present flow solver (INS3D code) has been validated by computing fundamental fluid 
dynamics problems such as channel flow, the flow over a backward-facing step, and flow 
over a circular cylinder. Three-dimensional cases include flow over an ogive cylinder, flow 
through a rectangular duct, wind tunnel inlet flow, cylinder-wall juncture flow, and flow 
through multiple posts mounted between two plates. These results have been reported 
in our earlier publications, and will not be repeated here. 1 3 ’ -8 30 The most striking 
application of the present code is the flow simulation in the power head portion of the 
SSME. 31,32 In the remaining part of this section, several important results obtained to 
date will be presented. 

Background of SSME flow analysis 

For future scientific and commercial applications, an upgrade of the SSME power head 
has been under way which will substantially increase the operating margin and the engine 
durability. To achieve this goal without increasing the weight and size of the existing 
engine, it became essential to understand the dynamics of the hot-gas flow in the engine 
power head. Because of the complexity of the geometry, an experimental approach is 
extremely difficult as well as time consuming. Computational simulation, therefore, offers 
an economical alternative to complement the experimental work in analyzing the current 
configuration, and to suggest new, improved design possibilities. In the past few years, 
major milestones have been established from this effort. In this report, highlights of our 
initial task are presented. 

In the SSME staged combustion cycle, the fuel is partially burned at very high pressure 
and relatively high temperature in the preburners. The resulting hot gas is used to run 
the turbine and is then routed to the main injector where, along with additional oxidizer, 
it is injected into the main combustion chamber. The Reynolds number of the primary 
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flow in the manifold is of order 10 6 per inch. Because of the high gas temperature, the 
Mach number is less than 0.12. The flow is turbulent and is practically incompressible. 

Figure 1 illustrates the current arrangement for the SSME power head components. Hot 
gas discharged from the gas turbine enters the annular turnaround duct (TAD), and ex- 
periences a 180° turn before it diffuses into the fuel bowl. This assembly is called the hot 
gas manifold (HGM). The gas flows into the main injector through three transfer ducts 
on the left-hand side of the power head (fuel preburner side) and enters into the region 
of the main injector posts. On the right-hand side of the power head (oxidizer preburner 
side), there are two transfer ducts connected to the right-hand side of the main injector 
assembly. Around the bottom portion of each liquid oxygen (LOX) post in the main in- 
jector assembly, there are a number of small holes through which the hot gas flows into 
the main combustion chamber. There it mixes with the oxidizer, which comes through 
circular passage along the centerline of each LOX post. As a part of the current redesign 
effort, a CFD study has been conducted to simulate the dynamics of the hot-gas flow in 
the porver head. 

The computer mod el and the grid 

A computational model of the power head is chosen to analyze critical areas where dynam- 
ics of the hot-gas flow is expected to have a significant effect on the overall performance 
of the SSME. As shown in Figure 2, the model starts from the gas turbine exit on the 
fuel preburner side, and extends to the main injector assembly. The main injector consists 
primarily of a bundle of LOX posts, which is physically modeled by a porous media. 

Figures 2a and 2b demonstrate the 3-D grid for the SSME HGM. They show a 
horizontal and a vertical cross section, respectively, of the HGM. Figure 2c illustrates the 
details of H-grids for the cross-section of the three transfer ducts. This H-grid is generated 
for a unit circle. Near the boundary the grid lines are concentric circles except in the 
vicinity of the four singular points. Using the nearly orthogonal grid in this unit circle, 
one can obtain H-grid for tubes or ducts of any given shape and dimension by a simple 
linear transformation. Figure 2d shows an unwrapped surface of the annular iuel bowl with 
openings. The elliptical transfer duct shown in Figure 2c represents an advanced two-duct 
design, which will be explained later in this section. 
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The grid for the entire HGM system is generated by using algebraic functions, and is written 
with a high degree of flexibility for changing geometric configurations. By specifying the 
shape, the dimension, and the desired number of transfer ducts, a grid for a variety of new 
HGM configurations can be obtained in a short time. The ducts described in this paper 
are connected directly to the fuel bowl without any fairings, while in the current engine the 
three transfer ducts are connected smoothly to the annular fuel bowl with fairings. This 
configuration with an abrupt change in geometry is more demanding computationally than 
smooth configurations. 

Multiple-zo ne computation 

A large number of mesh points is required to resolve the 3-D turbulent flow in the SSME. 
To overcome the limitation in computer core memory, the domain of interest is divided into 
several zones. This requires a special treatment at the interface for a smooth continuation 
of the solution between zones. Figure 3 illustrates a five zone arrangement for the HGM 
flow field. Zone 1 is allocated for the TAD and fuel bowl. Zones 2, 3, and 4 are for the 
three transfer ducts. The racetrack of the main injector is represented by Zone 5. Also 
shown in the figure are some overlapping grids in the various zonal interfaces. The grid 
is chosen to be continuous and smooth across the zonal boundaries. In this paper, the 
racetrack (Zone 5) is not included in the computation. Since the vertical plane through 
the center of the fuel bowl and the main injector is taken to be a plane of symmetry, only 
half of the HGM is computed. 

The equations of motion given in equations (2) and (7) are of hyperbolic type with 
parabolic-type viscous diffusion terms. Waves are propagating in both up- and down- 
stream directions while the solution approaches a steady state. The interfaces between 
zones of the present problem are locations where the geometry changes abruptly. There- 
fore, in the neighborhood of those interfaces, flow is expected to experience a rapid change. 
To maintain a smooth continuation of the solutions across these zones, and hence to achieve 
a stable and fast -converging computation, a means of providing adequate communication 
for the traveling waves must be established. Overlapping regions and a proper zonal in- 
terpolaton scheme are thus required for this purpose. 

A forward or backward differencing, if applied to the interfaces of multiple zones, would 
distort the geometric representation. To maintain a smooth transition of the flow field 
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across a zonal boundary, the Jacobian and the metrics at the interfaces are computed using 
grid points in neighboring zones. Then the pressure and the velocities, Q, are updated 
explicitly at each iteration. Let values at (n 1/2) denote the state of conditions to be 
used to advance the computation to n + 1. The values of Q n+1/2 for Zone 1 at the exit 
plane are obtained from the values of the corresponding plane of Zone 2 at n, i.e., 
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And values of Q" +1 / 2 for Zone 2 at the zonal interfaces are taken from the latest computed 
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When more than two points are overlapped, the latest values in the interior of this over- 
lapping region must be properly transmitted to the next zone. There are a number of ways 
to treat this problem. The simplest one is take an average of the two values computed in 
Zone 1 and 2, as below: 
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A scheme using updated zonal boundary values, but with no interior updates, has also 
been tested. Either way, converged steady-state solutions have been obtained. However, 
the scheme with interior updates converges at a much faster rate. 


Grid-induced error 


To generate grids for realistic 3-D problems like the SSME, it is convenient to combine 
different types of grids, depending on the problem’s geometric characteristics. In the 
present application, for continuity and smoothness across zonal boundaries, an H-type grid 
is chosen for the circular transfer duct as shown in Figure 2c. This is then is connected to 
the side-wall grid of the HGM as illustrated in Figure 2d. In generating this grid for the 
SSME, errors are introduced mainly due to grid singularities, skewness, and stretching. 

The present AF algorithm integrates the difference equations along the trans- 
formed coordinates, £, //, and ; directions. At the junction of the two H-grid directions, 
flow particles in the two coordinate directions could communicate only indirectly via the 
interior mesh points. This produces some corner-effect error. Even though this error is 
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not as severe as the one caused by external flows, an ad hoc method of eliminating this 
corner effect is devised based on a finite-element concept. Let j and k denote the indices 
for the grid points along the increasing ( and r/ directions, respectively. And let j = k = 1 
be the corner point, which is singular. First, the pressure at this point, pn, is determined 
by an extrapolation along the diagonal direction j = k. Second, P13 and P3] are obtained 
in a usual manner. Then p J2 and p 2 1 are established by an interpolation along the circular 
surface. 


The full viscous term given by equation 11a can be simplified to equation lib when the grid 
is orthogonal. Even though full viscous terms can be used, it is convenient and economical 
to keep only the orthogonal part. It is, therefore, of practical interest to estimate overall 
error caused by orthogonal formulation when a nonorthogonal grid system is used. Using 
the notation defined in Figure 4 a, the skewness can be related to the metric terms as 

dsi ■ ds t = (i ] x 2 + Vy 2 )[-j ) 2 


ds 2 ■ ds 2 = (t x 2 + t y 2 ){ d -j ) 2 

w d(dr t 
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As a quick measure of ail overall error caused by grid skewness, a 2 -D channel flow is com- 
puted using two grids, namely, 1) stretched Cartesian grid (orthogonal) and 2 ) nonorthog- 
onal grid where the skewness is controlled on the upper hall of the channel as sketched 
in Figure 4 b. In the computation only orthogonal terms are kept. Converged solutions 
on the lower half, where the two grids are identical, are then compared as shown in Fig. 
4 b. Total error depends additionally on the Reynolds number and the third-directional 
skewness. However, this quick experiment indicates that the orthogonal assumption can 
be used without significantly impacting on the overall solutions if the grid is reasonably 
orthogonal. 
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Turbulence models 


Several levels of turbulence models have been implemented into the code. These include 
the Baldwin and Lomax 33 algebraic model in which length scale is determined by the 
location of the maximum moment of vorticity. This model has been widely applied in 
the external flow problems. However, the maximum moment of vorticity is not as well 
defined for fully turbulent internal flows as for external flows. In particular, the moment 
of vorticity is almost constant for a fully developed pipe or channel flow except in the 
sublayer region. For the present problem, it is proposed that the length scale is determined 
by the point of minimum vorticity. This length scale is incorporated into an extended 
P ran d tl - K armann mixing-length theory. The combination of these automatically account 
for curvature effect. Full details of this model are given in other reports. 32,34 For the 
present SSME flow computations, this extended Prandtl-Karmann mixing-length model is 
used. A review on various levels of turbulence models caivbe found in Reference 35 or 36. 

Computed results 

In this part, flow solutions in a variety of different HGM configurations at various Reynolds 
numbers are presented. Here, the Reynolds number is based on the mean velocity and the 
duct width at the entrance of the TAD. 

Three-ci rc ular-duc t HGM 

First, steady-state solutions are obtained for the current three-circular-duct HGM at. 
Re- 1000. As illustrated in Figure 5, in the present analysis, the racetrack of the main 
injector is not connected to the ducts. The three transfer ducts are assumed to discharge 
the flow separately. There is no communication of the pressure between the center and 
the outer ducts at their exit planes. For this reason, small residual waves have remained 
in the computed results. However, the root mean square value of the change in the flow 
variables, AQ, has dropped below 10 3 , and an essentially steady-state solution has been 
obtained. 

In Figures 6a and 6b, velocity vectors are shown in t he horizont al and vertical cross-sections 
corresponding to Figures 2a and 2b. The flow in the center t ransfer duct, as illustrated in 
Figure 6b is highly nonuniform. and a large separation region is formed just downstream 
of the entrance to the transfer ducts. By comparison of the vector length in Figure 6a, 
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ihe flow in the center duct is much slower than it is in the outer ones. The results shown 
in these figures agree qualitatively with the airflow test data conducted at a Re of about 
10 6 . The predicted mass flow through the center duct is 9.8% of the total mass flow, which 
agrees with the test data. 

Figures 7a and 7b illustrate the 3-D velocity vectors at the unwrapped center plane and at 
a plane near the inner wall of the fuel bowl. A reverse-flow pattern is clearly visible near the 
inner wall. Three-dimensional swirl patterns are predicted in the vicinity of the entrance 
to the transfer ducts. Figure 7c is a photograph that indicates, by means of surface-streak 
(shear-pattern) visualization, the similar swirls at the corresponding locations in the airflow 
test . 


The existence of the swirls can be explained as follows: The flow coming from below has 
a large momentum due to the relatively small width of the annular duct. Among the 
streamlines of this flow in between the two ducts there exists a dividing streamline. This 
streamline has a stagnation point at the top of the fuel bowl as shown in Figures 7a and 
7b. On the left-hand side of this dividing streamline, flow is bent leftward to the center 
duct. Because of symmetry, a rightward flow is also approaching the center from the other 
side. When these opposite currents approach each other, another dividing streamline is 
formed with a stagnation point again at the top wall. The stagnation pressure forces the 
streams to bend downward, and at the same time, the streams make a right-angle turn 
into the circular duct. Conservation of momentum thus requires the formation of swirls. 

The pattern of the swirl and its center depends on the relative strength of the appoaching 
currents. Near the center duct, double swirls of equal strength are formed because of 
symmetry. In the vicinity of the entrance to the outer duct, the current approaching 
leftward from the rear part of the bowl is more massive than the one approaching rightward. 
A stronger swirl is thus formed and is located sideways toward the weaker stream. 

Figures 8 and 9 are the perpendicular cross-sectional views showing three different sect ions 
of the transfer ducts; namely, near the entrance, at the midsection, and near the exit 
plane. Near the entrance (Figures 8a and 9a), the velocity vectors in the center duct have 
symmetric double swirls, while the outer duct has a strong swirl accompanied by a much 
weaker one. The swirling velocities are largely reduced at the midsection and are physically 
dissipated before entering the main injector regions. 
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New two-elliptical-duct HGM 


From this computational flow analysis and also from experiments, the center duct of the 
current three-duct HGM is found to transfer a limited amount of mass flow (about 10% of 
the total flow). Also the transverse pressure gradient, remains large together with a large 
bubble of separation after the 180° turn. To improve the quality of the flow, a large- area, 
two-duct, design concept has been developed. In addition, the ducts are chosen to have an 
elliptical shape in order to distribute the mass flow evenly to the main injector region. An 
outer-surface grid for a two-duct model is illustrated in Figure 10. 

First, to reduce the large separation bubble after the 180° turn, a parametric study is 
performed to find the best possible configuration. In Figures 11a and lib, comparison 
of the current, three-duct configuration and the new two-duct design is illustrated. As 
shown in Figure 11a, a large separation bubble existing in the present design is practically 
removed in the new configuration shown in Figure lib. This is confirmed by experiment, as 
shown in Figures 11c and lid for the current and new designs, respectively, where velocity 
measurements at five different locations across the channel between the inner and the outer 
wall are shown. 

Figure 12a illustrates an example of the triple-swirl pattern in the elliptical duct near 
the entrance at Re=10 3 . Here, because of the absence of the center duct, the stream 
approaching from the left-hand side is much stronger than in the previous case. Another 
point of int erest is that the upcoming stream entering the elliptical duct direct ly from below 
is also more massive. The three currents are almost of the same strength, resulting in a 
triple-swirl flow. The swirling is greatly dissipated along the duct. Figure 12b shows the 
remaining small swirling vectors at the duct. exit . A steady-state, turbulent-flow solution 
for an HGM with two elliptical transfer ducts at Re - 10 5 has been obtained. 32 In Figure 
1 .‘ 5 , the swirling flow in the elliptical t ransfer duct is illustrated. In comparing these results 
with the laminar solution in Figure 12, it is seen that only a double-swirl pattern exists. 

The most significant aspect of the present study is to pinpoint the locations where flow 
experiences the most energy losses. An important measure of the energy losses is the 
mass- weighted average total pressure along the flow. Figure 14 illustrates the decreasing 
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coefficient of the mass- weighted total pressure along the centerline of the TAD, the fuel 
bowl, and the transfer duct. The total pressure coefficient C po is defined as 

n PO ” P03 


where 



P + + v 2 + ir 2 ) 

& 


dvr 


The discontinuities shown in the figure correspond to the entrance of the duct where energy 
fluxes are computed over different planes. In the figure, three different HGM configurations 
are compared. The initial two-duct design shows 28% less total pressure drop compared to 
the current three-duct version. After fine-tuning the two-duct configuration computation- 
ally, the pressure drop decreased even further to 36% less than the original configuration. 
This final configuration is then tested using cold air flow, which shows 40% reduction in 
pressure loss. 


CONCLUDING REMARKS 

This paper presents a summary of incompressible Navier-Stokes flow-solver development 
work. The SSME power head has been simulated using this flow solver. Computational 
results are favorably compared with test data, and offer information not readily available 
from experiments. The results show that CFD can reduce the development time and cost 
by suggesting the best possible configurations for final verification by experiments. For 
example, in redesigning the TAD, over 20 different configurations were studied computa- 
tionally, thus providing the best geometry to designers. Further study is in progress, and 
the total performance improvement will be compared in the future. At the time the present 
SSME analysis was performed, the comput ational model of the power head was designed to 
obtain solutions within a reasonable turnaround time. Therefore, the total number of grid 
points was limited and the model could include the TAD and transfer ducts only. Despite 
its limitation, the present application provides an excellent, example how the present CFD 
and computer capabilities can be integrated into the aerospace design process. 
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Figure 1. SSME power head component arrangement 
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Figure 2. (Concluded), c) H-grid for circular and elliptic cross section 
of transfer duct; d) unwrapped surface of annular fuel bowl 
(cross section C-C). 
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Figure 3. Multiple-zone arrangement 
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Figure 4. Grid effect on channel flow solution at Re=1000. 

a) Definition of grid skewness; b) Relative error due to skewness. 
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Figure 5. Inner and outer surface grid for a three-duct HGM. 
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(b) 


Figure 6. Computed velocity distribution at Re=1000. 

a) Top view; b) Vertical cross section of center transfer duct. 
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Figure 7. Velocity vectors on unwrapped surfaces. 

a) Unwrapped center plane; b) near the inner wall; 


c) experiment. 
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Figure 10. Outer-surface grid for a two-duct HGM. 


41 



cj 180 225 270 315 0 45 90 315 180 225 270 315 0 45 90 135 

> CIRCUMFERENTIAL LOCATION, 6 , deg CIRCUMFERENTIAL LOCATION, 0, deg 

(c) (d) 


Figure 11. Turn- A round- Duct redesign. 

a) Computed velocity vectors for the current design; 

b) Computed velocity vectors for the new design; 

c) Experimental velocity head measurements for the current design 
(6 defined in Figure 2a); 

d) Experimental velocity head measurements for the new design. 
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(b) 

Figure 12. Velocity vectors at cross sections of transfer duct in two-duct HGM (Re=10 3 ) 
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(b) 

Figure 13. Velocity vectors at cross sections of transfer duct in two-duct. HGM (Re=10 5 ) 
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Figure 14. Pressure losses in three-duct and two-duct HGM. 
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